编程讨论第一贴：讨论97年A题的编程方法

dcyu

A题 零件的参数设计
一件产品由若干零件组装而成，标志产品性能的某些参数取决于这些零件的参数。零件的参数包括标定值和容差两部分。进行成批生产时，标定值表示一批零件该参数的平均值。容差则给出了参数偏离其标定值的容许范围。若将零件参数视为随机变量，则标定值代表期望值，容差通常规定为均方差的3倍。
进行零件参数设计，就是要确定其标定值和容差，这时要考虑两方面因素：
一．当各零件组装成产品时，如果产品参数偏离预先设定的目标值，就会造成质量损失，偏离越大，损失越大；
二．零件容差的大小决定了其制造成本，容差设计的越小，成本越高；试通过如下的具体问题给出一般的零件参数设计方法。
粒子分离器某参数（记作y）由7个零件的参数（记作x）决定，经验公式为：
<img border=0 src=http://vip.6to23.com/dcyu/sxjm/shumo/gongshi.bmp onload="javascript:if(this.width>screen.width-300)this.width=screen.width-300">
        目标值（记作 ）为1.59。当y偏离  0.1时，产品为次品，质量损失为1000元； 当y偏离  0.3时，产品为废品，损失为9000元。
        零件参数的标定值有一定的容许变化范围；容差为A, B, C三个等级，用与标定值的相差值表示。A等为 1%, B等 5%, C等为 10%。7个零件参数标定值的容许范围及不同容差等级，零件的成本（元）如下表（符号/表示无此等级零件）：
                                       
        标定值容许范围        C等        B等        A等
        [0.075, 0.125]        /        25        /
        [0.225, 0.375]        20        50        /
        [0.075, 0.125]        20        50        200
        [0.075, 0.125]        50        100        500
        [1.125, 1.875]        50        /        /
        [12  , 20 ]        10        25        100
        [0.5625,0.935]        /        25        100
现进行成批生产，每批产量1000个，在原设计中7个零件参数的标定值为  容差均取最便宜的等级。
请你综合考虑y偏离 造成的损失和零件成本，重新知己零件参数，包括标定值和容差，与原设计比较，总费用降低了多少。


大家试试看做一下，用多种算法实现。

dcyu

没人愿意做一下，好象人气不够啊。

路过

好恐怖，看着就头痛，让我想一想，用条件遍历。

dcyu

偶贴一下自己的程序，网格法。
/*
  Filename : 97A.c

  Author dcyu

  Date 2002.9.13

*/

#include &lt;stdio.h&gt;
#include &lt;conio.h&gt;
#include &lt;math.h&gt;
#include &lt;stdlib.h&gt;

#define PI 3.141519265359

float f(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
return ((
      174.42*x1*pow(x3/(-x1 + x2),0.85)*
     sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7)))/x5
       );
}


float ff(float x1,float x2,float x3,float x4,float x5,float x6,float x7)
{
return((
      174.42*x1*pow(x3/(-x1 + x2),0.85)*
     sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7)))/x5
       );
}

float g1(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
        return 148.257*x1*x3*sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7))
                /pow((x2-x1),2)/pow(x3/(x2-x1),0.15)/x5 +
                174.42/x5*pow(x3/(x2-x1),0.85)*sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7));

}

float g2(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
        return -148.257*x1*x3*sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7))
                /pow((x2-x1),2)/pow(x3/(x2-x1),0.15)/x5 +
                87.21*x1*pow(x3/(x2-x1),0.85)*
                (0.792288*x4*pow( (1-0.36*pow(x4/x2,-0.56)) , 0.5)
                /pow(x2,2)/pow(x4/x2,0.4) + 3.0392*x4*pow(x4/x2,0.16)*
                pow( (1-0.36*pow(x4/x2,-0.56)) , 1.5)/pow(x2,2)        )/x5/x6/x7/
                sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7));

}

float g3(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
        return 148.257*x1*sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7))
                /(x2-x1)/pow(x3/(x2-x1),0.15)/x5 ;

}

float g4(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
        return -87.21*x1*pow(x3/(x2-x1),0.85)*
                (0.792288*pow( (1-0.36*pow(x4/x2,-0.56)) , 0.5)
                /x2/pow(x4/x2,0.4) + 3.0392*pow(x4/x2,0.16)*
                pow( (1-0.36*pow(x4/x2,-0.56)) , 1.5)/x2 )/(x5*x6*x7*
                sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7)));

}

float g5(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
        return ((
      -174.42*x1*pow(x3/(-x1 + x2),0.85)*
     sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7)))/pow(x5,2)
       );

}

float g6(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
        return -87.21*x1*pow(x3/(x2-x1),0.85)*
                (1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/x5/pow(x6,2)/x7/
                sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7));

}


float g7(float x1, float x2, float x3, float x4, float x5, float x6, float x7)
{
        return -87.21*x1*pow(x3/(x2-x1),0.85)*
                (1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/x5/x6/pow(x7,2)/
                sqrt((1 - 2.62*pow(1 - 0.36*pow(x2/x4,0.56),1.5)*
          pow(x4/x2,1.16))/(x6*x7));

}


float fai(float x,int l)
{
float rm,rn,rp,rt,pq,rrt,x1,x2,y,s,h;
float p1,p2,q1,q2;

if(x==0) return 0.5;
rp=1;
if(x*l&lt;0) rp=-1;
x1=fabs(x);
x2=x1*x1;
y=1/sqrt(2*PI)*exp(-0.5*x2);
rn=y/x1;
if(rp&gt;=0&&fabs(rn)&lt;1e-12)         return 1.0;
if(rp&lt;=0&&rn==0)   return 0.0;

rrt=3.5;
if(rp==-1) rrt=2.32;
if(x1&lt;=rrt)
{
         x1*=y;
         s=x1;
         rn=1.0;
         rt=0;
         do
         {
                 rn+=2;
                 rt=s;
                 x1*=x2/rn;
                 s+=x1;
         }
         while(s!=rt);
         pq=0.5+s;
         if(rp==-1) pq=0.5-s;
         return pq;
}
q1=x1;
p2=y*x1;
rn=1;
p1=y;
q2=x2+1;
if(rp!=-1)
{
         rm=1-p1/q1;
         s=rm;
         rt=1-p2/q2;
}
else
{
         rm=p1/q1;
         s=rm;
         rt=p2/q2;
}
while(1)
{
         rn++;
         if(rm==rt||s==rt) return rt;
         s=x*p2+rn*p1;
         p1=p2;         p2=s;
         s=x*q2+rn*q1;
         q1=q2;         q2=s;
         s=rm;
         rm=rt;
         rt=1-p2/q2;
         if(rp==-1) rt=p2/q2;
}

return rt;

}


char rongcha(int n)
{
        switch(n)
        {
        case 0: return &#39;A&#39;;
        case 1: return &#39;B&#39;;
        case 2: return &#39;C&#39;;
        default: exit(1);
        }

}

void main()
{
    FILE *fp;
        int i[7],cha[7],j,k,m;
        float s[3][7] = {{25, 50, 200, 500, 50, 100, 100}, {25, 50, 50, 100, 50, 25, 25}, {25, 20,
                      20, 50, 50, 10, 25}};
    float a[7] = {0.075, 0.225, 0.075, 0.075, 1.125, 12, 0.5625};
    float b[7] = {0.125, 0.375, 0.125, 0.125, 1.875, 20, 0.935};
    float tc[7] = {0.005, 0.03, 0.01, 0.01, 0.15, 1.6, 0.0375};
    float rong[3] = {0.01, 0.05, 0.1};
    float c[7],x[7],tx[7],ct[7],tcha[7];
    float y,mincost,w,t1,t2,cy,tcy,ty,cost,tst,lt,tlt,temp;

        mincost=100000000;
        m=5;
    if((fp=fopen(&quot;output.txt&quot;,&quot;w&quot;))==NULL)  exit(1);
        for(cha[0] = 1; cha[0] &lt;= 1; cha[0]++)
        for(cha[1] = 1; cha[1] &lt;= 2; cha[1]++)
        for(cha[2] = 0; cha[2] &lt;= 2; cha[2]++)
        for(cha[3] = 0; cha[3] &lt;= 2; cha[3]++)
        for(cha[4] = 2; cha[4] &lt;= 2; cha[4]++)
        for(cha[5] = 0; cha[5] &lt;= 2; cha[5]++)
        for(cha[6] = 0; cha[6] &lt;= 1; cha[6]++)
        for(i[0]   = 0; i[0]   &lt;= m; i[0]++  )
        for(i[1]   = 0; i[1]   &lt;= m; i[1]++  )
        for(i[2]   = 0; i[2]   &lt;= m; i[2]++  )
        for(i[3]   = 0; i[3]   &lt;= m; i[3]++  )
        for(i[4]   = 0; i[4]   &lt;= m; i[4]++  )
        for(i[5]   = 0; i[5]   &lt;= m; i[5]++  )
        for(i[6]   = 0; i[6]   &lt;= m; i[6]++  )
        {
                cost=0.0;
                for(j = 0; j &lt; 7; j++)
                {
                        x[j]=a[j]+(b[j]-a[j])*i[j]/m;
                        ct[j]=s[cha[j]][j];
                        cost=cost+ct[j];
                }
                c[0]=g1(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                c[1]=g2(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                c[2]=g3(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                c[3]=g4(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                c[4]=g5(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                c[5]=g6(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                c[6]=g7(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                cy=0.0;
                for(j = 0; j &lt; 7; j++)
                {
                        temp=c[j]*x[j]*rong[cha[j]];
                        cy=cy+temp*temp;
                }
                cy=0.333333333*sqrt(cy);
                y=f(x[0],x[1],x[2],x[3],x[4],x[5],x[6]);
                t1 = fai((1.6-y)/cy,1)-fai((1.4-y)/cy,1);
                t2 = fai((1.8-y)/cy,1)-fai((1.2-y)/cy,1);
                lt = (9000-1000*t1-8000*t2);
                w=1000*(cost+lt);
                if(w&lt;mincost)
                {
                        mincost = w; tcy=cy; tlt=lt;
                        for(j=0;j&lt;7;j++)
                        {
                           tcha[j]=cha[j];
                           tx[j]=x[j];
                        }
                        tcy = cy; tst = cost; ty = y;
                        printf(&quot;mincost: %f       y=%f\n&quot;,mincost, y);
                        fprintf(fp,&quot;mincost: %f       y=%f\n&quot;,mincost, y);


                }
        }
               
        printf(&quot;\n\nMin cost: %f , x:&quot;);
        fprintf(fp,&quot;\n\nMin cost: %f , x:&quot;);
        for(j=0;j&lt;7;j++)
        {
                printf(&quot;x[%d]=%f &quot;,j,x[j]);
                fprintf(fp,&quot;x[%d]=%f &quot;,j,x[j]);
        }
        printf(&quot;\nRongchadengji: &quot;);
        fprintf(fp,&quot;\nRongchadengji: &quot;);
        for(j=0;j&lt;7;j++)
        {
            printf(&quot;%c &quot;,rongcha(cha[j]));
                fprintf(fp,&quot;%c &quot;,rongcha(cha[j]));
        }
        printf(&quot;\ny = %f, cy = %f, cost = %f, lt = %f\n&quot;,ty, tcy, tst, tlt);
        fprintf(fp,&quot;\ny = %f, cy = %f, cost = %f, lt = %f\n&quot;,ty, tcy, tst, tlt);
        getch();
        getch();
               

}

路过

是呀，二天前开的论坛，人气高才怪呢。继续遍历...宣传...
源程序不错呀！版主一定是一个看整理资料的人

dcyu

其实这道题我认为是在全国竞赛中算是比较难的。关键是题意的理解。
而且编程也十分困难，对于算法来说，用上面的网格算法和Monte Carlo算法相对比较容易
编程，但是我看到有优秀论文用了模拟退火编程，退火算法的程序并不容易。
有人在训练的时候编过吗？

dcyu

人气不高没有关系，等以后人气旺的时候在来讨论也成。

路过

写错了，是[爱整理资料的人]。退火算法，用起来不爽。

dcyu

就是用google搜索“数学建模”排名第一的网站中有97a题的模拟退火算法，程序写得确实非常不错。我想很多人建模的过程不喜欢用退火的一个原因是退火算法中参数的确定很难，需要在运行过程中根据CPU和内存的实际情况制定参数方案。今年美国竞赛中就是自己不愿意用退火算法而采用了另外更弱的算法解才未能有好成绩的，所以编程的心得就是在赛前一定要很好对每一种算法进行实践，比赛中才能游刃有余。

Echo

有在86年B的紧急救护里用过退火法.
我们组的程序员研究了很久.
我看到十个小组中, 只有两个小组在用,包括我们这组.